home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / os2 / sysb091a.zip / sysbench / src / pmb_linpack.c < prev    next >
C/C++ Source or Header  |  1994-11-05  |  27KB  |  1,232 lines
  1. /*
  2. Translated to C by Bonnie Toy 5/88
  3.  
  4. You MUST specify one of -DSP   or -DDP     to compile correctly.
  5. You MUST specify one of -DROLL or -DUNROLL to compile correctly.
  6. You MUST specify a timer option(see below) to compile correctly.
  7.  
  8. To compile double precision version for Sun-4:
  9.    cc -DUNIX -DDP -DROLL -O4 clinpack.c
  10.  
  11. To compile single precision version for Sun-4:
  12.    cc -DUNIX -DSP -DROLL -O4 -fsingle -fsingle2 clinpack.c
  13.  
  14. To obtain   rolled source BLAS, add -DROLL   to the command lines.
  15. To obtain unrolled source BLAS, add -DUNROLL to the command lines.
  16.  
  17. PLEASE NOTE: You can also just 'uncomment' one of the options below.
  18. */
  19.  
  20. /* #define SP     */
  21. #define DP
  22. #define ROLL
  23. /* #define UNROLL */
  24.  
  25. /***************************************************************/
  26. /* Timer options. You MUST uncomment one of the options below  */
  27. /* or compile, for example, with the '-DUNIX' option.          */
  28. /***************************************************************/
  29. /* #define Amiga       */
  30. /* #define UNIX        */
  31. /* #define UNIX_Old    */
  32. /* #define VMS         */
  33. /* #define BORLAND_C   */
  34. /* #define MSC         */
  35. /* #define MAC         */
  36. /* #define IPSC        */
  37. /* #define FORTRAN_SEC */
  38. /* #define GTODay      */
  39. /* #define CTimer      */
  40. /* #define UXPM        */
  41.  
  42. #ifdef SP
  43. #define REAL float
  44. #define ZERO 0.0
  45. #define ONE  1.0
  46. #define PREC "Single "
  47. #endif
  48.  
  49. #ifdef DP
  50. #define REAL double
  51. #define ZERO 0.0e0
  52. #define ONE  1.0e0
  53. #define PREC "Double "
  54. #endif
  55.  
  56. #define NTIMES 100
  57.  
  58. #ifdef ROLL
  59. #define ROLLING "Rolled "
  60. #endif
  61.  
  62. #ifdef UNROLL
  63. #define ROLLING "Unrolled "
  64. #endif
  65.  
  66. //#include <stdio.h>
  67. #include <math.h>
  68.  
  69. static double st[8][6];
  70.  
  71. double pmb_linpack ()
  72. {
  73.    static REAL aa[200][200],a[200][201],b[200],x[200];
  74.    REAL cray,ops,total,norma,normx;
  75.    REAL resid,residn,eps;
  76.    REAL epslon(),kf;
  77.    double t1,tm,tm2,dtime();
  78.    static int ipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
  79.  
  80.    lda = 201;
  81.    ldaa = 200;
  82.    cray = .056; 
  83.    n = 100;
  84.  
  85. /*
  86.    fprintf(stdout,ROLLING);fprintf(stdout,PREC);
  87.    fprintf(stdout,"Precision Linpack\n\n");
  88. */   
  89.    
  90. //   fprintf(stderr,ROLLING);fprintf(stderr,PREC);
  91. //   fprintf(stderr,"Precision Linpack\n\n");
  92.  
  93.     ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
  94.  
  95.     matgen(a,lda,n,b,&norma);
  96.     t1 = dtime();
  97.     dgefa(a,lda,n,ipvt,&info);
  98.     st[0][0] = dtime() - t1;
  99.     
  100.     t1 = dtime();
  101.     dgesl(a,lda,n,ipvt,b,0);
  102.     st[1][0] = dtime() - t1;
  103.     total = st[0][0] + st[1][0];
  104.  
  105. /*     compute a residual to verify results.  */ 
  106.  
  107.     for (i = 0; i < n; i++) {
  108.            x[i] = b[i];
  109.    }
  110.     matgen(a,lda,n,b,&norma);
  111.     for (i = 0; i < n; i++) {
  112.            b[i] = -b[i];
  113.    }
  114.     dmxpy(n,b,n,lda,x,a);
  115.     resid = 0.0;
  116.     normx = 0.0;
  117.     for (i = 0; i < n; i++) {
  118.            resid = (resid > fabs((double)b[i])) 
  119.      ? resid : fabs((double)b[i]);
  120.            normx = (normx > fabs((double)x[i])) 
  121.      ? normx : fabs((double)x[i]);
  122.    }
  123.     eps = epslon((REAL)ONE);
  124.     residn = resid/( n*norma*normx*eps );
  125.    
  126. /*   printf("   norm. resid      resid           machep");
  127.    printf("         x[0]-1        x[n-1]-1\n");
  128.    printf("%8.1f      %16.8e%16.8e%16.8e%16.8e\n",
  129.       (double)residn, (double)resid, (double)eps, 
  130.            (double)x[0]-1, (double)x[n-1]-1);
  131.  
  132. fprintf(stderr," times are reported for matrices of order %5d\n",n);
  133. fprintf(stderr,"      dgefa      dgesl      total       kflops     unit");
  134. fprintf(stderr,"      ratio\n");
  135. */
  136.     st[2][0] = total;
  137.     st[3][0] = ops/(1.0e3*total);
  138.     st[4][0] = 2.0e3/st[3][0];
  139.     st[5][0] = total/cray;
  140.  
  141. //   fprintf(stderr," times for array with leading dimension of%5d\n",lda);
  142.    print_time(0);
  143.  
  144.     matgen(a,lda,n,b,&norma);
  145.     t1 = dtime();
  146.     dgefa(a,lda,n,ipvt,&info);
  147.     st[0][1] = dtime() - t1;
  148.     
  149.     t1 = dtime();
  150.     dgesl(a,lda,n,ipvt,b,0);
  151.     st[1][1] = dtime() - t1;
  152.     total = st[0][1] + st[1][1];
  153.     
  154.     st[2][1] = total;
  155.     st[3][1] = ops/(1.0e3*total);
  156.     st[4][1] = 2.0e3/st[3][1];
  157.     st[5][1] = total/cray;
  158.  
  159.     matgen(a,lda,n,b,&norma);
  160.     
  161.     t1 = dtime();
  162.     dgefa(a,lda,n,ipvt,&info);
  163.     st[0][2] = dtime() - t1;
  164.     
  165.     t1 = dtime();
  166.     dgesl(a,lda,n,ipvt,b,0);
  167.     st[1][2] = dtime() - t1;
  168.     
  169.     total = st[0][2] + st[1][2];
  170.     st[2][2] = total;
  171.     st[3][2] = ops/(1.0e3*total);
  172.     st[4][2] = 2.0e3/st[3][2];
  173.     st[5][2] = total/cray;
  174.  
  175.     ntimes = NTIMES;
  176.     tm2 = 0.0;
  177.     t1 = dtime();
  178.  
  179.    for (i = 0; i < ntimes; i++) {
  180.            tm = dtime();
  181.       matgen(a,lda,n,b,&norma);
  182.       tm2 = tm2 + dtime() - tm;
  183.       dgefa(a,lda,n,ipvt,&info);
  184.    }
  185.  
  186.     st[0][3] = (dtime() - t1 - tm2)/ntimes;
  187.     t1 = dtime();
  188.  
  189.    for (i = 0; i < ntimes; i++) {
  190.            dgesl(a,lda,n,ipvt,b,0);
  191.    }
  192.  
  193.     st[1][3] = (dtime() - t1)/ntimes;
  194.     total = st[0][3] + st[1][3];
  195.     st[2][3] = total;
  196.     st[3][3] = ops/(1.0e3*total);
  197.     st[4][3] = 2.0e3/st[3][3];
  198.     st[5][3] = total/cray;
  199.  
  200.    print_time(1);
  201.    print_time(2);
  202.    print_time(3);
  203.  
  204.     matgen(aa,ldaa,n,b,&norma);
  205.     t1 = dtime();
  206.     dgefa(aa,ldaa,n,ipvt,&info);
  207.     st[0][4] = dtime() - t1;
  208.     
  209.     t1 = dtime();
  210.     dgesl(aa,ldaa,n,ipvt,b,0);
  211.     st[1][4] = dtime() - t1;
  212.  
  213.     total = st[0][4] + st[1][4];
  214.     st[2][4] = total;
  215.     st[3][4] = ops/(1.0e3*total);
  216.     st[4][4] = 2.0e3/st[3][4];
  217.     st[5][4] = total/cray;
  218.  
  219.     matgen(aa,ldaa,n,b,&norma);
  220.     t1 = dtime();
  221.     dgefa(aa,ldaa,n,ipvt,&info);
  222.     st[0][5] = dtime() - t1;
  223.  
  224.     t1 = dtime();
  225.     dgesl(aa,ldaa,n,ipvt,b,0);
  226.     st[1][5] = dtime() - t1;
  227.  
  228.     total = st[0][5] + st[1][5];
  229.     st[2][5] = total;
  230.     st[3][5] = ops/(1.0e3*total);
  231.     st[4][5] = 2.0e3/st[3][5];
  232.     st[5][5] = total/cray;
  233.  
  234.    matgen(aa,ldaa,n,b,&norma);
  235.    t1 = dtime();
  236.    dgefa(aa,ldaa,n,ipvt,&info);
  237.    st[0][6] = dtime() - t1;
  238.  
  239.    t1 = dtime();
  240.    dgesl(aa,ldaa,n,ipvt,b,0);
  241.    st[1][6] = dtime() - t1;
  242.  
  243.    total = st[0][6] + st[1][6];
  244.    st[2][6] = total;
  245.    st[3][6] = ops/(1.0e3*total);
  246.    st[4][6] = 2.0e3/st[3][6];
  247.    st[5][6] = total/cray;
  248.  
  249.    ntimes = NTIMES;
  250.    tm2 = 0;
  251.    t1 = dtime();
  252.    for (i = 0; i < ntimes; i++) {
  253.       tm = dtime();
  254.       matgen(aa,ldaa,n,b,&norma);
  255.       tm2 = tm2 + dtime() - tm;
  256.       dgefa(aa,ldaa,n,ipvt,&info);
  257.    }
  258.    st[0][7] = (dtime() - t1 - tm2)/ntimes;
  259.    
  260.    t1 = dtime();
  261.    for (i = 0; i < ntimes; i++) {
  262.       dgesl(aa,ldaa,n,ipvt,b,0);
  263.    }
  264.    st[1][7] = (dtime() - t1)/ntimes;
  265.    total = st[0][7] + st[1][7];
  266.    st[2][7] = total;
  267.    st[3][7] = ops/(1.0e3*total);
  268.    st[4][7] = 2.0e3/st[3][7];
  269.    st[5][7] = total/cray;
  270.  
  271.    /* the following code sequence implements the semantics of
  272.       the Fortran intrinsics "nint(min(st[3][3],st[3][7]))"   */
  273. /*
  274.    kf = (st[3][3] < st[3][7]) ? st[3][3] : st[3][7];
  275.    kf = (kf > ZERO) ? (kf + .5) : (kf - .5);
  276.    if (fabs((double)kf) < ONE) 
  277.       kflops = 0;
  278.    else {
  279.       kflops = floor(fabs((double)kf));
  280.       if (kf < ZERO) kflops = -kflops;
  281.    }
  282. */
  283.    if ( st[3][3] < ZERO ) st[3][3] = ZERO;
  284.    if ( st[3][7] < ZERO ) st[3][7] = ZERO;
  285.    kf = st[3][3];
  286.    if ( st[3][7] < st[3][3] ) kf = st[3][7];
  287.    kflops = (int)(kf + 0.5);
  288.  
  289. //   fprintf(stderr," times for array with leading dimension of%4d\n",ldaa);
  290.    print_time(4);
  291.    print_time(5);
  292.    print_time(6);
  293.    print_time(7);
  294. //   fprintf(stderr,ROLLING);fprintf(stderr,PREC);
  295. //   fprintf(stderr," Precision %5d Kflops ; %d Reps \n",kflops,NTIMES);
  296.    return kflops;
  297. }
  298.      
  299. /*----------------------*/ 
  300. static print_time (row)
  301. int row;
  302. {
  303. /*fprintf(stderr,"%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n",
  304.        (double)st[0][row], (double)st[1][row], (double)st[2][row], 
  305.        (double)st[3][row], (double)st[4][row], (double)st[5][row]);
  306. */
  307. }
  308.       
  309. /*----------------------*/ 
  310. static matgen(a,lda,n,b,norma)
  311. REAL a[],b[],*norma;
  312. int lda, n;
  313.  
  314. /* We would like to declare a[][lda], but c does not allow it.  In this
  315. function, references to a[i][j] are written a[lda*i+j].  */
  316.  
  317. {
  318.    int init, i, j;
  319.  
  320.    init = 1325;
  321.    *norma = 0.0;
  322.    for (j = 0; j < n; j++) {
  323.       for (i = 0; i < n; i++) {
  324.      init = 3125*init % 65536;
  325.      a[lda*j+i] = (init - 32768.0)/16384.0;
  326.      *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
  327.       }
  328.    }
  329.    for (i = 0; i < n; i++) {
  330.       b[i] = 0.0;
  331.    }
  332.    for (j = 0; j < n; j++) {
  333.       for (i = 0; i < n; i++) {
  334.      b[i] = b[i] + a[lda*j+i];
  335.       }
  336.    }
  337. }
  338.  
  339. /*----------------------*/ 
  340. static dgefa(a,lda,n,ipvt,info)
  341. REAL a[];
  342. int lda,n,ipvt[],*info;
  343.  
  344. /* We would like to declare a[][lda], but c does not allow it.  In this
  345. function, references to a[i][j] are written a[lda*i+j].  
  346. */
  347.  
  348. /*
  349.      dgefa factors a double precision matrix by gaussia